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ia LNT ROBUCTION 


The parabolic partial differential equation 


= = y-(a(x,t,u) Vu) + £(x,t,u, Vu) (x,t) © 2x(0,t] 
has been the subject of the development of several numerical 
approximation techniques. The equation is known as the 
heat equation in a specialized form and is useful in the 
study of heat transfer and thermodynamics in general. It 
also has applications in other fields of engineering and 
science, such as fluid dynamics and meteorology. 

The development of Galerkin methods for the solution 
of a general parabolic and hyperbolic equation has been 
relatively recent [3,6,9]. This thesis describes a program 
for the solution of the general parabolic problem using the 
central different Laplace modified alternating direction 
Galerkin scheme described in [5,6,8]. The program was built 
in steps beginning with the one dimensional case and then 
extended to the two dimensional case. 

The program solves the problem for a rectangular domain 
Q, and uses nine Gauss quadrature points for integration in 
the plane. It also allows for nonlinearities in the forcing 
PiimetwiOMmeri%,t,Uu,V7u) and the function a(x,t,u) can be 
extended to have nonlinear terms of Vu. The theoretical 
PeeOLtcworesStaoll1ty and COnvergence are available for these 


cases in [4]. 





The program was written in the Fortran IV computer 


language and tested on an IBM 360/67 computer system. 


Me DeEetNi TIONS, NOTATION, AND PRELIMINARY CONCEPTS 


In the following section several concepts and definitions 
which are essential for the understanding of the thesis are 
presented. Two specific spaces of functions will be consid- 
ered and definitions of the norms in the spaces will be 
given. 

Let (eke with smooth boundary 02, be the domain of 
interest for the function spaces. The first of these spaces 
is that of square integrable or measurable functions, L*(Q). 


Lt weL?(Q), then the norm of w is defined as 
Iw | =f w ax < 
ele (2) ° 
L* (Q) is a Hilbert space with respect to the inner product 
Sec si(x) g(x) dx 
Q 
The Sobolev space H° (Q) 1s equivalent to RE RON) 
The second space considered will be the Sobolev space 


H+ (9), defined as 


{w|DPw eb (a), Ip| = 0,1} 


2 
I 


where 


pP = —___*_______ , Ip (p. ) 


| 
ws 


10 , 





The norm of w is defined as 


7 Z a dW 2 We 
Ill dca) = CMwll aca) * 8 Mosel paca) 


Of special interest for the development of this paper 


is the space Ho a subspace of He defined as 


Ho = {w|w eH" (Q) aad sw |= 0, £or all x €909)} 


The norm for this subspace is defined to be 


Nn 
OW 2 Ly 
IW ledcgy = £ 2 Weer p2 (9)! 
Hy (2) gay OXy 0 LEED) 


A major concept which 1S used in the formulation of 
Galerkin's method is that of the weak solution to a differ- 
ential equation. Let L be the differential operator 


describing the differential equation, that is 
Lu = f (2re1) 


A function u is said to satisfy the weak form of equation 


(2.1), 1£ for every test function v in some test space 


Se vee ES (2-2) 


If u is not sufficiently differentiable to be substituted 
iMmCcwi jut does Satisfy (2.2), then u is said to be a 


weak solution. 


Ie P 





To illustrate the weak solution, consider 


Lu = es UL. = F.(scet.) , pee 2 = (0,2) tt > 0. 
mee eye — tt ),ct) = O° F u(x,0) = Up (x) ; Sere 1 
(2453) 
The weak form of this equation is 
2 = i 

Sue War 7 = <f,v> fi @ieare 11 v € H(t) (2.4) 

mene) = uc l,t) = 0 

u(x,0) = Up (x) 
Integrating by parts results in 

a a 
Seo Meg ct. Suey os 2B ieee FOr va ee v € Ho (Tt) 
tt, 0) = Up (x) (25) 


Note that u € Ho (I) canmoecee Solution to (2.5), that 
1s a weak solution of (2.3), however a strong solution of 


(2.3) must satisfy ue C(I) for eacn t. 


EZ 





Ili. THE FORMULATION OF THE GALERKIN METHOD 


Rieownormuitatien OL Galerkin methods for the solution of 
differential equations is based on the weak form of the 


equation. Consider the nonlinear parabolic equation 


oF 


— i eaoateinevimiet £(xX,t,u,Va) , (x,t) 6 2x(0,T) 


u(x,0) = Up (x) Pee, Se (SecA) 


eae) 6 00x {(0,T] 


| 
© 


an, Cc) 


and the test space Ho (2) . Ciearly for any Vv Ho (2) 


= 


ou 


= an V°(a(x,t,u) Vu)v + £(x,t,u,Vu)v 


Taking the integral over the region & yields 


ii iy Geel (alx,t,u) Vu)v dx + J £(x,t,u,Vu)v dx 
Q Q Q 


Using Green's first formula for integration by parts in 


multiple dimensions results in 


i! sty gx = / a(x,t,u) Vu v dx - f a(x,t,u)VurVv dx 
Q Q Q 


Pemex tu. vu)y ax ©. 
Q2 


13 





Suace vy = 0 on of it follows that 


(i eicioawk wujsvu Vvedx = ..0 
022 


and 


i ~ v ax + f a(x,t,u) Vu°Vv dx = A fester, VL) Ny oa x 
Q Q 


Therefore u satisfies 


—a(x;0),v> = Sy v7 fOr val) VE Hy (2) 


Let M be a finite dimensional subspace of Hg (2) . The 
continuous-time Galerkin method is to find for each te (0,T] 


a differentiable map U(-,t):{0,T] +M , such that 


<P /V> + <a(x,t,U)VU,VV> = <£(x,t,U,VU) ,V> (3.3) 


Vv > = <Up 1 V> , t=0, Vem 


Now let M be a subspace of Ho (2) such that 


ue 


M = Span (Ww) +Wo7---/W,) , where {w i=1 


fi is a linearly 


independent set. Under these conditions it is possible to 


write 


fms 


U(x,t) = os (t) wa (x) ‘ 


j=l 
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Substituting in equation (3.3) and letting V = Ws 


< 


lms 
Ios 


1, CE) wi (x) sw (x) > eee (SO, te, U ) 


; ig Mei) Ware (eg 1» Vw; (x) > 


J 1 


= <f£(x,t,U,VU) ,w; (x) > 


and 


< 


Shay SUN ed = <Up -W; (x) > 


lis 


al 


Since integration is a linear operation and the inner product 


defined is not over the time space, it follows that, 


n n 
ag (t) <wes (3) wy (2x) > + Y a. (t)<a(x,t,U) Vw. (x) , Vw, (x) > 
= <f£(x,t,U, VU) ,w, (x) > (3.4a) 
and 
n 
om eg A) als) 7 = <Ug 1W,; (x) > il, 2 eee OC eb) 


Equations (3.4a) and (3.4b) can be written as a system 


of n-nonlinear ordinary differential equations 
Maceo a att) = F 


3 o5)) 
Ma(0Q) = b 


Is 





where S(a) = (S55) and M = (m5 5) are nth order matrices and 


si4 SELL Bo Ell) Ne ea) Ae) 


Il 


Ms 4 Sue (x) ,W. (x) > 


and 


F = Ci.) ? ft. aa <f(x,t,U,VU) ,w,; (x) > 


1 
Ail 
Qa = (G5 1Go ree er) 
pa oe we >, <u. , Wa Se Ww =o 
Gee Qt 2 7! OF 


Using the properties of the functions War it can be 
shown that the matrix M is positive definite, and with the 
additional property that a(x,t,u) is bounded above and below, 
the matrix S can also be shown to be positive definite. If 
in addition a(x,t,u) iGreen yi) oat isty Lipscn1 eZ 
conditions, it can be shown that the system of differential 
equations has a unique solution. 

In the one dimensional case, an approximate solution to 
equation (3.5a) can be obtained by combining it with a 


finite difference scheme to discretize in time as follows. 


Pesumang, Lor simplicity, that a(x,t,u) = a = constant 
and f(x,t,u,u,) = £(x, ¢) 
igs tot m+1l ) om ay 
M(————-) + S(=—,——) = F''® (3.6) 


is 





This is a system of n equations and n unknowns and can be 
Esmplified to 


1 
(ameats)ya™t? = (am-ats)a™ + 2atr™? 


Or 


[Ope 


1 
(2M+AtS) (a me 2Aesa + 2Acr 


erwin e 


m m+% 
° =~ e s = < 
where 0. So Mie and F £ (t+ 3 


In the above equation the time variable t has been discretized 
cO oe = mAt , where t = T/N , N a positive integer. 

This method is known as the Crank-Nicholson-Galerkin method. 
The one dimensional program was written to set up the 
matrices defined above and solve the algebraic equation (3.7). 

A minor alteration in the program would allow the use of 
(3.7a), and give the advantage of reducing the computational 
roundoff error. The program is also limited to the linearity 
Simplifications made above. 

The required accuracy of the method in the two dimensional 
case leads to a large number of basis functions for the space 
and consequently would require excessively large matrices. 

The matrices have a band structure and the band width of the 
Matrices is dependent on the number of intervals used in the 
grid set up for the region. These large matrices can be 

avoided by implementing an alternating direction method, and 


using a tensor product basis for the space. This tensor 


Jy) 





product basis reduces the problem to working with matrices 
Mie! Nave a Constant band width of seven, in the case of 


M being Hermite cubic splines. It also greatly simplifies 


the problem when the region 2 R? 


1s rectangular. 
Let Le = eg SF be a basis for 
one dimension of the space M and Jy = Span(n, (y) ,n5(y),.-- 


a (y)) be a basis for the other dimensions of M. Me and 


M, are bases for one dimensional subspaces of Ho (1,) and 
a _ = 
Ho (ty). such that M = gu = oi igi Uy eal We 


eee ne ) . Also define 


x Y 


<£,9>, = f fg dx 


I 
Xx 


Sie Ge edi 
Ys 5 
y 


Since Mae forms a basis for the space M, we can write 


Nn 1g 
x 


y 
U(x,y,t) = . a hae Ve (6, 0) Q Ng fy)) 


There are several time discrete methods with which to 


proceed at this point. The program has been designed to use 


the alternating direction version of the centered difference 


Laplace modified Galerkin method described in [5,6,8] and 


defined as 


18 





m+] au 


U ge 





m m 
= oe pve <a (x,t,U) VU , W> 
+ 2“ at< a” 57um av, = <f™ (x GHUy vue v > 
aX OY r 9X OY evr Ve f 
and 
<U,V> = <Ug1V> V eM, a M, bcd) 


where the superscript of the functions denotes the time 


step at which it is to be evaluated, and 


y2ym . ym l ~ 2y™ 4 yt 
and 

2 

omens : L? (Q) 


oXd0Y 


A is subject to the restriction 


1 elle 
Dose F 4max = 7 Max alee (5-1 ) ) 
x 2 
t>0 


mem Stability [5,6]. 
Equation (3.7) can be written in the algebraic tensor 


pereoauct form 


m+1 


C+ 2AAtA_) @ (C_ +2XAtA = 2(2AA QA +A AC 
(C., x) BIC, y) a (2A (C, MA, +A MC) 


2 2 m 
+ 4)° (At) (A, @AL)) a 
i 


2 2 m~ 
ae (Cac, 2hAt(C, @A, +A, BC) 4“ (At) (A, @AL)) a 


+ 2arom 


v9 7 





where 


m m m m 
oe ee eesytey) VU V(E., ee tos (x,t,U,VU),&,, se 
and 
2 liter a 1 ee 
ag SEG 2 Ay = SUE 
Define e” = ofl giles , then equation (3.7) can be 


written as 


(C,#2AMtA,) @ (Cy#2AAtAL) (eT -e™) = -2(C,@C,)e™ + 200" 
(3.8a) 
and 
Ce (ere ar oe (3. 8b) 


The scheme (3.8) adds significant efficiency to the program 
and also reduces roundoff error. 

The Laplace centered difference method initially requires 
two sets of the coefficients, those at time zero and at one 
time step later. At time zero the set of coefficients are 
cbtained by an L (2) Peoyecte rons oO: thes anitial conditions 
into the space, but the second set presents a problem where 


a variety of methods are available. Two of these methods 


20 





were attempted, expanding the initial conditions in a 
Taylor's Series about t = 0 and using the Laplace forward 
difference method. The former involved taking the deriva- 
tives of the initial condition function and presented a 
problem when these derivatives resulted in zero. The method 
which was finally used was the Laplace forward difference 
scheme. It is defined as follows 


ae m 


a) = (At) om 


(C,. 1 AtA,) @ (Cc, a AAtA,,) (a (34:9) 


At first glance the matrices in (3.8a) and (3.9) seem to be 
different, however in the case of this method A is under 


the restriction 


Max a7, Wl) So Crm ade aU 


>’ 
Vv 
NR] 
ey) 
I 
NS) 


thus a prudent choice of A results in identical matrices and 
it is not necessary to form a second set of coefficient 
Matrices in the program. It is also noteworthy that 6" has 
the same form as in the previous method and the subroutine 
used to compute o™ may be used for both with a correction of 
> for the forward difference scheme. These advantages ied 


me che final choice. 


The tensor product of matrices has the property that 


(B@D) = (B QI) (I @D) 


Zl 





This property enables a significant simplification of two 
aspects of the programming problem. First, the tensor 
product multiplication on the right hand side of the equation 
can be accomplished by computing first 


(IQ Cyjen = y" 


and then 


m 

(C, a iy 

By properly ordering the vector e” this results in a series 
of multiplication by the relatively small matrices Cy eigyel (Gas 


Y 
that is, matrices which have a constant band width of seven 


hor Sesame YE cubics rather than a band width which depends 

on the number of intervals in the grid of the problem as in 

C,ac,. 
This same idea can be used in soiving the linear system. 

Defining the right hand side of the equations (3.8,3.9) 

to be ay let 


pede m m 


(IQ (C,, + 2AAtA,)) (e e) = Y 


Using this, solve first 
hail 
((C, + 2Adta,) Ci, = 6 


and then solve 


(1 @ (Cc, + 2AMtA,)) (es Ses ci 


Ze : 





Again proper ordering of the vectors involved results in 
the simplification of solving a system of one dimensional 
problems rather than the larger problem of solving the full 
Pemnsor product matrix. 

The two dimensional program has been designed to solve 
the problem in the manner described above. The basis 
functions used in the program are Hermite bicubic splines, 
taken as the tensor products of one dimensional Hermite 
cubic splines, and will be described in Section IV. MThe 
procedure described above has the drawback that it is in 
general good only on rectangular polygons, and must be 
altered significantly to handle regions other than simple 
rectangles. It can be changed to handle regions which are 


unions of rectangles Sale 


92 





ieee EB CUBIC SPLINES 


The Galerkin procedure is dependent on a set of linearly 
independent basis functions for the space considered in the 
problem. Hermite bicubic splines provide such a basis. 
These splines are a high powered tool for the approximation 
Se smooth functions. 

The Hermite bicubic splines are taken as the tensor 
product of one dimensional Hermite cubics. These cubics 
have a piecewise polynomial nature, where the polynomial is 
of degree less than or equal to three. Hermite cubics have 
globally continuous first derivatives. Thus, the tensor 
megcuct Of the one dimensional cubics has the property that 
both first partials and the first mixed partial derivatives 
are continuous. These properties are sufficient for the 
approximations desired as a result of the program [1]. 

The package of subroutines described in [2] was used for 
the program. It consists of subroutines which compute the 
values of the B-splines at specific points, and given a set 
of coefficients, computes the value of the approximated 
function at specified points. The high degree of flexibility 
of the package was the major reason for its use in the 
program. The program can be extended to use higher degree 
B-splines which gives a greater degree of smoothness and 
Bceuracy tO the approximating function. Also, the extension 


of the program to higher dimensions can be accomplished by 


24 





taking a higher degree tensor product, and making the 
appropriate alterations in the computational procedures. 


Gonsiaer the interval [0,1], and the partition of the 


interval A: ¢ — Xo Ss xX = pia ‘s ae =]. cer h; = KS es ee 
meet. = [X. -,xX.]. Define h = max h. , and h = min h. 
1 a- i" 1 3 af = ‘i aL 


The Hermite cubics are piecewise cubics over r. 130) gamers yt Oa gen iar 
and are members of ane HiVewOrder Or accuracy for Hermite 


cubics 1S given by the following well-known theorem: 


Theorem. If ne cme On on, ee een 


~ 4—- 
inf llu-all| siz) < Ch” Mull pags) 6 OS 8S 3 
ue 


where M is the space of Hermite cubics. 


25 - 





ve SekeROR Teo TIMATES 


The usefulness of a procedure to obtain an approximate 
solution to a differential equation iS meaSured by its 
degree of accuracy and the efficiency and costs of its 
implementation. The purpose of this section is to derive 
an error bound for the continuous time Galerkin method for 
a linear version of equation (3.1). The following theorem 
ema its proof is found in [8]. It is restricted to the case 
Pere a(x,t) is not a function of u. Proofs for the nonlinear 


case are also available [3,6,9]. 


Theorem 1. Let u be the solution to 
ou : _ 1 
<aE ,v> + <a(x,t)Vu,Vv> = 0 Vv € Ho (2) ee eee Ore | 
(5a) 
<u(x,0)> = <Up (x) > xen 
and let U be the solution to 
oU 7 _ 
<aa ,V>+<a(x,t)VU,VV> = 0 ect weet © GO, ]| (a2) 
<U,V> = <ug,/V> t = 0 
and 
O< Co < aay (xr t) < cy (5 430 


26 





Then there exists constants C and dg, which depend on 


Cor Cys ane. ep such that 


2 ic > 
o<tct Il u-Ull5 2 (0) (eta : Herold 2 ke 
max ~1) 2 . ety 
£ Clo cee Ju-tl] £2 eo) (e) + J || u-a || nd (a) (t) dt 
+ a : —U 2 a CS 4 
; |] Se Curt) |] F279) (t) at) (5.4) 


where u is an arbitrary map of (0,T]J into M. 


Smeeere Replace v with V in (5.1) and subtract (5.2) from 


wo. 1) i. 


< 2 pie Saxe) Va VV =: 10 
dU = 
i wi POV os <ats, t) VU,VV> = 0 
results in 
< oien lew Gece, viu-Ue,yV> =—= 0 4 


Let e= u-U and let V = eté , where e= u-u , then 


< ,ete> + <a(x,t)Ve,Vete> = 0 


de 
ot 


27 





Note 


< se ,e> = f oe e dx 
Q2 
= 5 oF A a dx 
= 3 ae llell £29) 
and 
<aVe,Ve> > cp |lel| nL (9) from (5.3) 
and 
|<ave,Vé>| < |laVel] 2g) * IlVéll 29) 


[A 


Se, |lell) 2 seme ed 
1 ell whee) + W&ll yLeo 


from Schwarz's inequality. 


Young's inequality states ||f|| - llgll < ¢ ll£ll* + = Ilgll* 


Using this inequality we get 


2 Sele hey 2 
| <aVe, VE> | < eT Nell ut (a) a 4e | € |] Hg (2) 


Thus 


ia Z 9e ~ 2 
> at Ile || L2 (9) cs SSE ,e> + Co ell dca) 


p ear ay, 
see, Mlellit(a) * Te NEN a2 (a) te (0,T} 


Za 





and 
lle] | 2. 5 i oe e> + (c -€Ec ) lel| 2 
1a eat §)) ay oe 0 i) HA (2) 


R32 
res HEN F:2 (9) te (0,T} 


Zz 
Let ¢€= Cy) /2cy a ieleveio tee@te 16 = C) /C. 


d 2 0@ ~ 2 
ae | e|| ii) 2 Spe 02 ae Se Hell gd (a) 
Z 
Pacegrating from 0 to T 
ell 2 (t) - |le||72,,,(0) + ¢ 5, el] 71 dt 
L* (2) Lé (2) 05 Ha (2) 


18 


15 
~ || 2 0e@ ~ 
2 ee We dt -2 f <<= ,@> dt 
eee 0 Hy (82) 0 dt 
and 
e E iE = 
uf gues die aes - ff ey 9e, dat 


= 2 ~1| 2 
|<e,é>| < Elle |] 72 (9) + C Ell 5209) 
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iy <e,-—| dt < (s Heal a5, (JG) neche as a re : (GeySchs), 
at a L4 (2) eee es) 
me rollows that 
2 _ 2 : 2 
eV tte Hell £2 gy (0) 1 fo A Nell 2 (a) dt 


2 2 anes 
See iig2 cy St? * 8 fell ;2¢p) tO) + ; Hell ad cay oe 
. 3 9 ce 
. 
eee 
7 A sell y2(¢Q) at) 


tiws for ¢ sufficiently small 
2 ae - 2 5 
Hell £2 (9) t) om le || Hd (2) Ls 
< C cf || e 2 dt + |le| 2 (0) + re é - Gite 
a Y Well 22 (9)  L2(9) J MEM ad co) 
ey, wee ee 
+ HET] 52 (9) (7) 1 HET] £2 CQ (0) + 4 sell £2 ¢9) dt) 
Srenwall's inequality: Let £,g, and h be piecewise continuous 
bemmegative functions defined on am anterval a«< t <b and 
assume that for each t ¢« [a,b] 


T 
fie) + h(t) <“g(t) + fe>t(s) ds 
_ 0 


30 





Then 


ie 
f(t) + h(t) < s eSg(s) ds + g(t) 
a 
and if g is nondecreasing 
eae 
=e) fF h(t} ve g(t) 


Thus 
2 : 2 
| e|| ee (Q) (t) ai 0 i | e]| Ha (2) dt. 
max Z Saas 
< © (gegen ell 52 (9) (©) ~ : Nell ga (a) at 


—_ 


Aen ey, 
Now by the definition of U(x,0) 
|| ¥-U || L (9) (0) a || u- | L (9) (0) 


and the desired result is obtained, that is 


2 2 
pew =U) eee Ca) eee 4 e-Oil d ca) ot 
max eee e Ae 
C octet Ju-F]] 529) (et) + J eral 2 (ay dt 


t 
+ 2. (v= |] £2 (6) (e) at) 


Syl 





This implies that the error bounds are dependent on how well 
u and its first partiai derivative with respect to time can 
be approximated by the arbitrary map u. For Hermite cubics, 
) u-u || 2 is O(h*) , and |lu-dl| .1 is O(h?) 
L? (2) Ho (2) 

With further restrictions it can also be shown that the 
a (2) error is onan 

In [3] J. E. Dendy, Jr. proves the following theorem 


which gives error estimates for the discrete time case. 


mneorem 2. Let u be the solution of 


<4 Povo etn aka) VM, VVee— <E(*,uU) ,Vv> ; v € Hy (2) 
u(x,0) = up (x) ea ii) 
u(+,t) Ho(Q*> , te (0,7) 


and assume that A > i Cc CO Stee. 2166 16,0) 5 Soe We ell 


Ala a 
mee Ef = pty , where UV eM, and M satisfies 
i p-s 
inf le-Oll 53 (9) eon Hull pr eg) [3]. Let W be the 
mapping from: (0,T] >M defined by 
<a(-,u (u-W),V> + <u-W,V> = 0 for all VeM 
Ba gu d¢u 2 rc ' 
If Mere ee 6) Xx (OF) ); and ere 2) lin (Ch A Sa oI (0) es a Dae 


then for 1<N <M and h and At sufficiently small 


eZ 





2 2 N- 
9 E ffae™||* at + at ENT + et fi 


— 
me 


2 
+ 2d7at || oy (a. malo s eo ie) 
| wy) | | = SOE 2 S02 ocNAt hr , 
if 
2 
re one 3 1 2 4 
IE" Iy + WE Wy + 4t Igggy -BD)ITS < e(4t) 


For the case of Hermite cubics, r = 4, and the conditicns 
on M are satisfied. U is the approximate solution obtained 
in the Laplace modified alternating direction method. 

Using the program it can be shown that the order of the 
error is actually this by approximating the value of w where 
the error is O(h”). This approximation can be obtained in 
the following way: ~ 
Let ey be the error for h = hye and e. be the error for 


h= h then 


2 f 


Wd 


ae he 


eat oo 


2 
log(e,/e,) = Ww log(h,/h,) 


w = (log(e,) - log(e,))/(1og(h,) - log(h,)) 


oe 





Using this method and choosing t = oes we obtain the 
results presented in Table l. Table 2 shows similar results 
mor the time error. 

To obtain the estimates in the two tables the following 
problem was used for easily obtained analytic solution. 


Be Yay f Sg = 2(e ©) ((x-x?) ie (y-y7) ate 5 (x-x*) (y-y*) } 


u(-,t) = 0 fOiwall x cof 


u(x,y,0) = (x-x?) (y-y) 


In this test A= Z was used in the discrete time 
formulation. 
TABLE 1 
ea ye L-2 NORM MAX NORM 
Ore2 0100) 0.0400 0.0003680 0.0006815 
0.1000 0.0100 0.0000286 0.0000619 


THE VALUE OF W FOR THE MAX NORM IS 3.46006. 
Pie oben OF W FOR THE L-2 NORM IS 3.68510. 


TABLE 2 
HL DT L-2 NORM MAX NORM 
0.2000 0.0400 0.0003680 0.0006815 
0.1000 0.0100 0 .0000286 0.0000619 


THE VALUE OF W FOR THE MAX NORM IS 1.73003. 
tite Or Ww FOR THE L-2 NORM IS 1.84255. 
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Vi owe ST -PROELEMS 


The utility of the program is evident in several fields 
of engineering. The following test problems will illustrate 
its use in two of these. The problems are taken from the 


fields of fluid dynamics and heat transfer. 


Exe RAYLEIGH'S PROBLEM FOR A CORNER 
In the first problem we consider the diffusion equation 


in two dimensions with the following boundary and initial 


conditions 
is (1/Re) (uy, + uly) , ey 2.0 
u(x,y,0) = 1 
m(Or 7, tc) = u(x,0,t) = 0 


Wiaeeyerte) te cSt sc vail) yes co 


where Re is the Reynold's number of ae ie Wbinirel bietelias a 
consideration. 

This is known as Rayleigh's problem for a corner, and 
the solution describes the impulsive motion of a right 
angled corner formed by two infinite flat plates and is 
used to infer the steady flow along the corner, with leading 


edge at t=0. 
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Bote enc upuapose OL Numerical approximations, using a 
Reynold’s number of 1000, x=4 and y=4 can be considered 
Bomintinity and the problem can be stated for use in the 


program as 


ig ee ae 
u(x,y,0) = 1 

HGOr yee) = atx, O,t) = 0 
Mier) u(x,4,t) = 1 


and 


ave, U) =~ 171000 


£(x,y,t,U/U,/U,) = 0 
A nonuniform grid with a fine mesh near the x=0, y=0 
corner is useful for an accurate description of the boundary 
layer behavior near the corner. Illustration (1) is a 
diagram of the grid used. 
The exact solution to this problem can be found analyt- 


ically to be 
Wien, tc) = err (Xx) veri (y) 


where 
1/2 


ms 
i 


(x/2) (Re/t) 


yf) ieee 


< 
ll 


S16 





Ti 
Ril 


Beas io ee ee SS SER 


at “i 


al] 


Me peledeqsct hone.) 





The nature of the solution is that of a shock wave 
traveling through the rectangle. The error in the approxi~ 
mation seems to be greatest on the wave and small even when 
the grid size is quite large. The error was also larger at 
early times near the boundary discontinuity than later in 
the process. Table 3 is a sample of the output of the 
program for this problem. The output points were taken near 
the discontinuity on the boundary since it is this area 
which is of greatest interest. 

For this problem i’ = .00025, At = .03 , and the 


largest Ax was .75. 


B. THE HEAT TRANSFER PROBLEM 
The equation for heat conductance for heterogeneous 
isotropic solids and frictionless incompressible fluids is 


oc me = Vo (kK u) + Gg" 


where the solution u(x,y,t) represents the temperature of 


the material. op is the density of the material, c is the 
heat capacitance, k is the thermoconductivity, and q,)'" is 
an internal generation term. pc is a constant and can be 


combined with k and dg to give 


ou = V-(k' u) + q'"/pc 


linmciicmtorm k*° 1s the thermal diffusitivity. 
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For most materials the thermal diffusitivity can be 
approximated by a linear function of the temperature. It 


ms typically of the form 


oan = Kee a(u- Up) 
mamere K is the initial value of the thermal diffusitivity 
and a iS a small real number. 

The internal generation term is highly variable and 
dependent on the type of problem. A typical function 


encountered is of the form 


o8 (u/ug) 


~~ 


Analytic solutions for problems such as this are not 
available and it is in these cases that approximation 
techniques are invaluable. 

The problem for the program run will be as follows 
(u/ug) 


u, = V- ((100 + .005 (u-uy) ) Vu) +e 


u(x,y,0) = sin(mx) + sin(tTy) 
2 

MO, tc) = Ul, y7e)) = ee! C Sin(tTy) 
2 

aise.) OPS ieee ee GT Eee) aaaiees 


neve eae = LOO O05 a — (can (ix) + san (ny) )} 


)+si 
E(xsy, teu, Ur) = _(u/(sin(1x)+sin(Ty))) 
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with uniform grid, with Ax = Ay = he= 1/10 , time 
step At = .01 and XK = 30. 

Plots 1-8 show the general nature of the solution. 
Plot 1 shows the solution at an early time while it retains 
much of its initial form. The area near the boundary is 
beginning to become smaller. The internal generation term 
causes the oscillatory nature of the solution. After a 
longer period of time the solution begins to take on a 


steady state of zero. 
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Peon (7) “TIME =-eil 








= 2.55 


(8) TIME 


PLOT 








VII. CONCLUSIONS AND RECOMMENDATIONS 


Two major conclusions can be reached from the applica- 
tion of the program to the various problems. First the 
program iS in general more efficient than the standard 
implementation of Galerkin methods which achieve the same 
degree of accuracy. It is also possible to achieve a higher 
degree of continuity and approximation in the solution using 
higher degree B-splines. It is possible that changes not 
evident to the author may be made to improve the efficiency 
of the program. Most of the computational effort in a given 
time step iS in evaluating u’ and vu" at the gauss points; 
more efficient evaluations or placement of the quadrature 
points can possibly increase the speed of the program. It 
should be noted that the matrices [see (3.8),(3.9)] arising 
need be factored only once since they are time-independent. 
Second, optimal error beunds for the method have not been 
derived in the case that the Laplace modified forward 
difference method is used to obtain the second set of 
coefficients so that the centered difference method may be 
used, but this does not seem to have any effect on the 
Feuution in application. 

Using the program as a basis the Galerkin methods may 
be extended to a third dimension, using extension of the 
milcermiating Girection £6 obtaim it. The basic ideas used 
in the development of the program can be utilized to program 


the second order hyperkolic problem using the Galerkin methods. 
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APPENDIX A 


tEE COMPUTER PROGRAM 


A. THE PROGRAM VARIABLES 
Pepnam = the coefficients Of the Splines at time m. 
CX,CY,AX,AY - the matrices Cy Cyr Ay Ay respectively 
described in Section III 
meres the grid points on the x and y axes 
TKX,TKY - the knots system on the x and y axes, used 
for the de Boor routines 
GX,GY - the Gauss weights computed in WEIGHT 
EKCX,ETAY —- the values of the basis functions at the 
Quadrature points 
ALXR,ALXL,ALYB,ALYT - the coefficients for the projection 
of the boundary conditions at time m 
GEx7EPY — the gquadrature points on the x and y axes 
NX,NY - the ane: of intervals on the x and y axes 
TAX, I AY - the dimension of the matrices defined above 
Pie At 
TIME - time m 
B. THE MAIN PROGRAM 
The main program calls various subroutines to accomplish 
the necessary procedures. Brief indications of the function 


of each subroutine will follow. The program begins by 


setting up the grid desired and its corresponding knots 


eecem and then computes the necessary matrices. It stores 


ope 





parts of the matrices which will be changed but needed for 
future corrections. The first set of coefficients is 
obtained by projecting the initial conditions into the 
Space. The second set is computed using a forward differ- 
ence Laplace modified step. The program then steps forward 
in time using the central difference Laplace modified 


method. 


ee) MATCOM 

MATCOM computes the matrices for the program and stores 
them in band symmetric storage used by IMSL library routines. 
Zt also stores the quadrature points and the values of the 


basis functions at the quadrature points. 


Dee FMATCO 


FMATCO projects the initial conditions into the space. 


i. WEIGHT 
WEIGHT computes and stores the product of the quadrature 


wemehts and Ax and Ay for the appropriate interval. 


F, GRID 


GRID computes the grid points. 


Se. KNOTS 


KNOTS computes the knot system for the de Boor routines. 


H. UEVAL 
UEVAL computes the value of the function at specific 


points (because of inefficiency its use should be lime ted). 
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i MULTIP 


Mm 


MULTIP computes (CX CY) e in manner described in 


Section III. 


wee OOLVER 
SOLVER solves system of the form (CX I)(I_ CY) e™ = B 


in the manner described in Section III. 


K. ANEVAL 
PREVAL evaluates the function and its first partial 


derivatives at the quadrature points. 


life RTHDSD 


RTHDSD computes the right hand side of equation (3.8). 


ie ALAYBC 


ALXYBC projects the boundary conditions into the space. 


ieee f£UNCTION F 
Function F defines the necessary functions for the 


problem being solved. 


See JHE DE BOOR SUBROUTINES 


The de Boor subroutines perform various operations to 


obtain B-splines [2]. 


5 3 





APPENDIX B 


THE ee PROJECTION OF FUNCTIONS 


The te projection of the initial conditions and other 
functions involved in the program were accomplished using 
Gauss quadrature. Appendix A describes several subroutines 
in which this is done. To achieve the desired accuracy it 
Mase necessary to use nine quadrature points for each interval. 
In the subroutine WEIGHT, the Gauss weights are combined with 
the value of Ax and Ay for their respective intervals and 
Saved. This is necessary in the case of non-uniform grids. 
The quadrature points are stored for future use in the 
subroutine MATCOM. 

To project a function into the space, we must take its 
inner produce with each of the basis functions. This means 
we must know the values of the basis functions at the 
quadrature points. These are computed using the de Boor 
routines and saved in MATCOM. In the case of one dimensional 
B-splines there are only four non-zero basis functions on 
each interval and each basis function iS non-zero on no more 
than two intervals. Thus on each interval the quadrature is 
taken with the functions and each of the four non-zero basis 


-NMCtLOnNS. 


th Cha 


For example, on the I x-interval and the J y-interval, 
we evaluate the function at the nine quadrature points and 


take its product with the value of the basis functions at 
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the point and the appropriate Gauss weights, and then sum 
over the nine points in the interval for each of the four 
basis functions. This procedure is continued by moving to 
the next interval, and using the appropriate basis functions 


mor that interval. 
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LIST OF PROGRAM 
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